Computing Laplace Transforms for Numerical Inversion Via Continued Fractions

It is often possible to effectively calculate probability density functions (pdf’s) and cumulative distribution functions (cdf’s ) by numerically inverting Laplace transforms. However, to do so it is necessary to compute the Laplace transform values. Unfortunately, convenient explicit expressions for required transforms are often unavailable for component pdf’s in a probability model. In that event, we show that it is sometimes possible to find continuedfraction representations for required Laplace transforms that can serve as a basis for computing the transform values needed in the inversion algorithm. This property is very likely to prevail for completely monotone pdf’s, because their Laplace transforms have special continued fractions called S fractions, which have desirable convergence properties. We illustrate the approach by considering applications to compute first-passage-time cdf’s in birth-and-death processes and various cdf’s with non-exponential tails, which can be used to model service-time cdf’s in queueing models. Included among these cdf’s is the Pareto cdf. Keywords—computational probability, numerical transform inversion, continued fractions, Laplace transforms, S fractions, complete monotonicity, Padé approximants, cumulative distribution function, birth-and-death process, Pareto distribution. Many probability density functions (pdf’s) and cumulative distribution functions (cdf’s) of interest in queueing models and other probability models arising in operations research can be effectively computed by numerically inverting Laplace transforms; see Abate, Choudhury and Whitt [1], Abate and Whitt [4], [5] and references therein. The biggest challenge in this approach, when there is a challenge, is usually computing the required Laplace transform values, because convenient closed-form expressions for Laplace transforms often are not available. In this paper we point out that continued fractions can sometimes serve as a basis for effectively computing the required Laplace transform values needed in the inversion algorithms. A simple motivating example is the steady-state waitingtime pdf in the M/G/1 queue. The classical PollaczekKhintchine (transform) formula gives the Laplace transform of the steady-state waiting-time pdf in terms of the Laplace transform of the service-time pdf. Thus we can compute the waiting-time transform values in order to compute the waiting-time pdf or cdf by numerical inversion whenever we can compute the service-time transform values. A possible difficulty, however, is that we might want to consider service-time pdf’s for which convenient explicit expressions for the Laplace transform are unavailable. Indeed, this difficulty often arises when we consider distribuSubject classifications: Probability distributions: calculation by transform inversion. Queues, algorithms: Laplace transform inversion. mathematics, functions: Laplace transforms and continued fractions tions which have non-exponential tails, e.g., which cannot be represented as phase-type distributions. The present paper provides a way to address this problem: Under favorable circumstances, we may be able to construct a continuedfraction representation of the service-time Laplace transform that enables us to compute the service-time Laplace transform values, which in turn enables us to compute the waiting-time Laplace transform values needed to perform the desired numerical inversion. A specific example covered by this approach is the Pareto pdf. For background on continued fractions and their use for numerical computation, see Baker and Graves-Morris [12], Bender and Orszag [13], Chapter 12 of Henrici [26], Jones and Thron [28], Section 5.2 of Press, Flannery, Teukolsky and Vetterling [32] andWall [35]. Applications of continued fractions in statistics and applied probability are described in Bowman and Shenton [15] and Bordes and Roehner [14]. More recently, Guillemin and Pinchon [20], [21], [22], [23] have used continued fractions to analytically derive important properties of queueing models. A summary of that work is contained in Dupuis and Guillemin [16]. However, continued fractions evidently have not been suggested previously as a way to numerically compute transform values in order to perform numerical transform inversion. The use of continued fractions is an alternative to computation of Laplace transforms via infinite-series representations, which we recently discussed in Abate and Whitt [9]. We make an explicit numerical comparison to show that continued fractions can be far superior in some circumstances, even when the series converges geometrically. (See Section 6.) Here is how the rest of this paper is organized: In Section 1 we briefly define continued fractions and specify the basic recursive algorithm for numerical computation. In Section 2 we discuss the relation between continued fractions and power series. There we show how to compute the continued fraction elements from the moments of a probability distribution (which are related to the coefficients of a power series — the moment generating function). In Section 3 we point out that completely monotone pdf’s can be identified with special continued fractions called S fractions, which have nice convergence properties. In Section 4 we show how continued fractions can be used to compute the Laplace transforms of first-passage-times pdf’s in birthand-death processes. We can exploit S fractions for this purpose because first passage times to neighboring states have completely monotone pdf’s. The rest of the paper is devoted to numerical examples. In Section 5 we consider the M/M/∞ busy period, which is a special case of a first passage time in a birth-and-death process. In Section 6 we consider the beta mixture of exponential (BME) pdf’s from Abate and Whitt [8] [9] and show that continued fractions can be much more effective for computing Laplace transform values than the previously considered infinite-series representations. In Section 7 we show how the continued fractions associated with the BME pdf’s can be used to compute the Laplace transforms of other pdf’s related to the BME pdf’s, including a Pareto pdf. Finally, in Section 8 we discuss continued fraction representations of Laplace transforms of other pdf’s. 1. Continued Fractions An (infinite) continued fraction (CF) associated with a sequence {an : n ≥ 1} of partial numerators and a sequence {bn : n ≥ 1} of partial denominators, which are complex numbers with an = 0 for all n, often called elements, is the sequence {wn : n ≥ 1}, where wn = t1 ◦ t2 ◦ . . . ◦ tn(0), n ≥ 1 , (1.1) and tk(u) = ak bk + u , k ≥ 1 ; (1.2) i.e., wn is the n-fold composition of the mappings tk(u) in (1.2) applied to 0, called the n approximant. If limn→∞ wn = w, then the CF is said to be (properly) convergent and the limit w is called the value of the CF. We write w = Φn=1 an bn or w = a1 b1+ a2 b2+ a3 b3+ · · · . (1.3) When we consider Laplace transforms of pdf’s, the CF elements an and bn will be functions of the complex variable s. In particular, we will consider special CFs called S fractions (S for Stieltjes), which can be expressed as w ≡ w(s) = 1 1+ a2s 1+ a3s 1+ a4s 1+ · · · (1.4) where ak is real and positive for all k. However, S fractions may appear in other forms, because CFs have many equivalent representations. Indeed, for any sequence of complex numbers {cn : n ≥ 0} with c0 = 1, the CF Φn=1 cn−1cnan cnbn (1.5) has the same approximants as the CF in (1.3); see p. 478 of Henrici [26]. We call such CFs equivalent and use the notation . It is significant that there is a relatively simple recursion for calculating the successive approximants of a CF, due to Euler in 1737. In particular, given the CF in (1.3), wn = Pn Qn , (1.6) where P0 = 0, P1 = a1, Q0 = 1, Q1 = b1 and Pn = bnPn−1 + anPn−2 Qn = bnQn−1 + anQn−2 (1.7) for n ≥ 2. In performing numerical calculations, it is prudent to renormalize after, say, every 10 iterations by dividing the current values of Pk, Qk, Pk−1 and Qk−1 all by Qk, and then proceed. We will be interested in the special case of S fractions, as in (1.4). Based on our computational experience, we conclude that the S fraction converges rapidly and is easy to compute if an = O(1) as n → ∞. If an = O(n) as n → ∞, then the S fraction converges, but more slowly and requires more work to calculate. If the elements grow much faster, then computation is likely to be infeasible. (We will give examples later.) The S fraction with an = n k can be shown to be convergent if and only if k ≤ 2; see p. 486 of Henrici [26]. 2. Power Series and Continued Fractions Continued fractions are intimately related to power series. This relationship is useful in probability applications because the moments of probability distributions can be regarded as coefficients of a power series, namely, themoment generating function (mgf). Let f be a pdf on the nonnegative real line with associated cdf F (t) ≡ ∫ t 0 f(u)du, t ≥ 0, and associated complementary cdf (ccdf) F (t) ≡ 1 − F (t), t ≥ 0. Assume that the pdf f has finite moments of all orders, i.e.,